Levy area logistic expansion and 
simulation 

By Simon J. A. Malham and Anke Wiese 

Maxwell Institute for Mathematical Sciences and School of Mathematical and 
Computer Sciences, Heriot-Watt University, Edinburgh EHI4 4^S, UK 
O ' (30th June 2011) 

(N ■ 

We present a new representation for the Levy chordal area for a two-dimensional 
Wiener process conditioned on its endpoints. This is based on an infinite weighted 
sum of Logistic random variables. We develop a numerical simulation algorithm for 
the Levy area based on truncating this series and simulating the tail by a suitable 
Normal random variable. We show how to improve the efficiency of the algorithm by 
Cr . I approximating higher order terms, which are large sums of independent identically 

PLh ' distributed Logistic random variables, by two separate methods: using a suitable 

Normal approximation and, sampling directly from the fixed density function for 
the logarithm of the product of decimal magnitudes of independent identically 
distributed Exponential random variables. To implement a strong Milstein numer- 
ical integrator for a stochastic differential equation driven by a multi-dimensional 
Wiener process, we must maintain a local mean-square error of order the cube of 
the stepsize. To achieve this prescribed accuracy, the latter two Levy area sampling 
^ \ methods we propose, reduce the number of uniform random variables required to be 

sampled over each timestep, a measure of the complexity, from reciprocal square- 
root complexity to logarithmic complexity. 
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1. Introduction 



Fundamental and crucial to high order strong simulation of stochastic differential 

equations driven by a multi-dimensional Wiener process, is the need to accurately 

k><( \ and efficiently simulate the Levy chordal areas when the diffusion vector fields 

5_j ■ do not commute. Thusfar this has represented a substantial technical difficulty. 

Cd I In the leading implemented methods, for example Levy's Normal expansion with 

a suitable Normal approximation of the tail sum, the number of random variables 
required to accurately simulate the Levy chordal area in the Milstein method, scales 
like the inverse square-root of the stepsize (Wiktorsson 2001). Here we introduce 
three new simulation methods for the Levy area based on a new expansion in 
terms of Logistic random variables. Beyond a given accuracy threshold, all three 
new simulation methods deliver highly accurate Levy area samples with orders of 
magnitude less computational effort. Indeed for two of the methods, the accuracy 
achievable scales logarithmically with the computational effort. 

Consider the problem of strongly simulating a stochastic differential equations 
driven by two independent Wiener processes W} and W"^ . Assume that the cor- 
responding diffusion vector fields do not commute. To implement an order one 
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Figure 1. Sample Levy chordal area — this is the area between the two-dimensional 
Wiener process shown and the chord joining the endpoints. 

Milstein method we must, on each computational timestep of size h, generate two 
independent N(0, v/i) sample Wiener increments, AW^{h) and AW^{h), and also 
a sample of the Levy area defined by 

1 rt+h pTi 

^W:=2/ / d<dW^2 -dW^^^d<. 

By Stokes' Theorem this is the area enclosed by the two-dimensional Wiener path 
and the chord joining the endpoints of the path at t and t + h; see Figure [T] 

In his seminal 1951 paper. Levy introduced an equivalent definition of the 
chordal area as follows. Using the Fourier series representation for the Brownian 
bridges based on the Wiener processes W^ and W^ given AW^{h) and AW^{h), 
Levy showed that the chordal area has the series expansion 

^ k=l 

where C/fc, Vfe, X^^Yu are independent standard Normal random variables. More gen- 
erally orthogonal eigenfunction series expansions for Wiener processes are referred 
to as Karhunen-Loeve decompositions after independent work by Karhunen (1947) 
and Loeve (1945). In the same 1951 paper. Levy proved that the characteristic 
function (j) — (j){C) corresponding to the probability density function cj) for the Levy 
area A{h), given AW'^{h) and AW'^{h), is 

'^(0 = -|||^exp(~ia2(i/i^coth(i/.0 - 1)), 

where a^ = {{AW^{h)f + {AW^{h)f)/h. 

Kloeden, Platen & Wright (1992) suggested truncating the Levy's series ex- 
pansion above to include N terms, and using that to generate approximate sam- 
ples. The mean-square error associated such a truncation scales like h? /N . How- 
ever in an important advance, Wiktorsson (2001) proposed simulating the tail sum 

by a matched Normal random variable {h/^/2'n:) (tt^/G — X]fc=i V*^^) ^ where 
Z ^ N(0, 1). With this correction, the mean-square error scales like h"^ /N"^ . In fact 
Wiktorsson achieves much more than this, deriving an effective sampling method 
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for the Levy areas involved in a Milstein approximation of a stochastic differential 
equation driven by any high-dimensional Wiener process; see Gilsing & Shardlow 
(2007) for a practical implementation. In this paper we will restrict ourselves to a 
two-dimensional Wiener process. 

We measure the effectiveness of any such simulation methods as follows. To 
successfully implement the Milstein method, we must ensure that the mean-square 
local truncation error is of order h^ . Hence over each computational timestep, if 
we used Wiktorsson's method above, we would need to take N ^ hr^l"^ . We will 
measure the computational effort associated with such an approximation, by the 
number of uniform random variables e that are required to generate such a Levy 
area sample. We will call this the complexity of the simulation method; also see 
Ryden & Wiktorsson (2001). Roughly, for Wiktorsson's method, the number of 
uniform samples we require scales like N. Note, we ignore: (1) the negligible effort 
associated with converting the uniform samples to standard Normal samples; (2) 
constant multiplicative factors (there are four standard Normal samples required 
at each order) and (3) the negligible effort associated with sampling the tail. Hence 
for Wiktorsson's method, the complexity required to achieve accuracy of order h^ 
is e~^''^. The smaller the complexity, the more effective is the simulation method. 

Gaines & Lyons (1994) proposed Marsaglia's rectangle- wedge-tail method for 
generating Levy area samples. However, Wiktorsson (2001) observes that "this 
method is complicated to implement, however, and occasionally requires numer- 
ical inversion of the characteristic function". Gaines & Lyons followed this in 1997 
with an approximate simulation method based on replacing the Levy area by its 
conditional expectation on intervening Brownian path information (also see Lord, 
Malham & Wiese 2008 and Malham & Wiese 2008). The complexity of this latter 
method is £~^. 

Ryden & Wiktorsson (2001) developed a method based on an expansion for the 
Levy area in Laplace random variables. Observing that the characteristic function 
(j) is the product of the characteristic function of a Logistic random variable and a 
Poisson mix of Laplace random variables, they derived the representation 

Pk 






27rV ^^ k 

^ k=l 3 = 1 

where X ~ Logistic(l), and for each fc € N we have: Pk ^ Poisson(a^) and for 
j — I,. . . ,Pfe, Yjk ^ Laplace(l); aU independent. Note that the expectation of Pk 
is a^, independent of k. In a practical simulation we would truncate this series 
to include all the terms ioT k ^ N and simulate the tail sum by the Normal 
random variable (ft,/27r) a yJ^jkZ where Z ^ N(0, 1). This simulation method has 
complexity e^^'^. 

Some related recent papers of interest are as follows. In a promising paper 
Stump & Hill (2005) derived a series expansion for the Levy area density function 
based on analytically computing the inverse Fourier transform of the characteristic 
function 0. The terms in the resulting series expansion decay like 1/fc^ with k as 
the summation index (see page 407 in Stump & Hill). However they only provide 
comparisons of the distribution functions and do not implement a sampling method. 
Levin & Wildon (2008) have developed methods for generating moments of the Levy 
area utilizing combinatorial shufSc product structures. Also, assuming Hormander's 
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ellipticity condition, Cruzeiro, Malliavin & Thalmaier (2004) have shown how to 
orthogonally evolve the frame bundle associated with the underlying system to 
avoid computing the Levy area entirely. 

Let us now summarize the ideas underlying the simulation methods we propose. 
Our starting point is Ryden & Wiktorsson's expansion for the characteristic func- 
tion (j) of the Levy area. As already noted, the first factor in the representation for 
(f> above, is the characteristic function associated with a Logistic random variable. 
We thus focus on the second factor, and for the moment, the expression in the 
exponent. Using an elementary hyperbolic function identity which permits indefi- 
nite iteration, we derive an infinite series expansion for the exponent with terms of 
the form (^a^ 2") (z/2")/sinh(z/2") where n is the summation index. Applying the 
exponential function to this series for the exponent generates an infinite product 
of exponentials of such terms. Now we first observe that characteristic functions 
of exponential form correspond to a Poisson mix of random variables whose char- 
acteristic function is the exponent. Then second we observe that the exponents 
(the individual terms from the series) are characteristic functions corresponding 
to scaled Logistic random variables. Hence we can represent the Levy area, given 
AW'^{h) and AW'^{h), as the following series of Logistic random variables (see 
Theorem 12.11 below) : 






2tt\ ^—' 2'^ 

\ ri=0 fc=l 

where for n = 0, 1, 2, . . ., the P„ ^ Poisson(ia^2") are independent Poisson random 
variables, and for n — 0,1,2,... and k — 1,2,...,P„, X,Xn^h ^ Logistic(l) are 
independent identically distributed standard Logistic random variables (mean zero, 
variance 7r^/3). 

This representation is the basis of the three simulation methods we propose. The 
first simulation method is a simple truncation of this series representation together 
with a tail approximation. Note that at each order n, we must on average generate 
an increasing number, namely ia^2" uniform random variables. This means that the 
complexity of this method is e~^ without, and e"^'^ with, the tail approximation. 
Hence, in principle this does not represent an improvement in complexity over 
the Levy or Ryden- Wiktorsson methods, though in practice the computational 
effort required to achieve a given accuracy is far superior to the Ryden- Wiktorsson 
method, and beyond a high accuracy threshold, also superior to the Levy method. 
See Figure [H where we provide an explicit accuracy versus computational effort 
comparison. The second method we propose, the Normal approximation, recognises 
that at each order we must sum an increasing number of Logistic random variables. 
Invoking the Central Limit Theorem, we replace sums of Logistic random variables 
over a threshold number by the appropriate Normal approximation. This of course 
has a dramatic effect on the computational effort required, without in practice 
seemingly compromising the accuracy of the simulation method — see Figure [5] The 
third method we propose, the Exponential product approximation, recognises that 
a sum of P independent identically distributed standard Logistic random variables 
can be represented by the difference of the logarithms of products of P independent 
standard Exponential random variables. We derive an explicit asymptotic form for 
the density function for large P for the logarithm of the product of P Exponential 
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Levy and Logistic approximations 




- Levy 
Logistic 

- Logistic norm. 

- Logistic exp. 

- Levy w. taii 

- Logistic w. taii 

- Logistic norm. w. taii 

- Logistic exp. w. taii 



Figure 2. The absolute value of the difference between the sample variance and the true 
variance (1 + a^) h? /12 versus the CPU time required to compute the simulation with 
5 X 10* independent samples. Here we fixed ft = 1, so with random Wiener increments 
AW^(h) and AW^{h) for each sample, we set a^ = ¥.{{AW^(h)f + {/\W^{h)f} = 2. 

random variables, using the method of steepest descents. We utilize this explicit 
form to tabulate the inverse distribution function for P — 10^, 10^, 10^ and 10'"'. 
Combining samples from these tabulated inverse distribution functions with small 
sums (< 10^) of Logistic random variables generates the highly effective Exponential 
product approximation shown in Figure [2] 

Our paper is structured as follows. We derive our Logistic expansion for the 
characteristic function for the Levy area in §2 and provide an algorithm for its 
implementation including a tail approximation. Then in §3 and §4 we respectively 
derive the Normal and Exponential product approximations based on the Logistic 
expansion. In both sections we discuss the complexity of these two methods. Lastly 
in §5 we provide further simulation details and conclude. 



2. Logistic expansion 

We propose a different expansion to that considered by Ryden and Wiktorsson 
(2001) producing a series representation for A{h) with quite different properties. 

Theorem 2.1. The Levy area A{h) conditioned on the given Wiener increments 
/^W^{h) anrf AM^^(/i) is equivalent in distribution to the following series of Logistic 
random variables: 

h / oo P„ ^ 



27r 



n=0 



fc=l 
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where for n = 0,1,2, .. ., the Pn ^ Poisson{-^a^2^^ are independent Poisson random 
variables, and for n = 0,1,2,... and k = 1,2,...,P„, X,Xn^k ^ Logistic{l) are 
independent identically distributed Logistic random variables. Any X ^ Logistic{l) 
is efficiently simulated by X ^ log(C//(l — [/)) where U ^ Unif{[0, 1]). 

Proof. We begin by observing the following elementary identity for hyperbolic func- 
tions, cothz = cothz/2 — 1/sinhz. Iterated N times, this generates the identity 

N 



cothz = cothz/2^+1 - V . / ,„ 



n=0 



Hence we see that 

exp(— ia^(zcothz — 1)) 



n=0 

where we set 



= exp(ia^$:2" . -i/^).exp(-la^(.eoth(./2-+^) 
n exp(^ia22" . _i-^ jexp(-ia22") • exp(f^.(z, a)). 



En{z, a) := - \a^ (z coth (z/2^+i) - 2^+i) . 



Note that £n{z, a) ^- as A^ — ^ cx), for all a; € R. We have thus established that 
for all ^ e R, setting z = £,h/2, we have 

'^^('^'^^^ = ^h^eV^ • n ^ UinheV2"+i j ' 

where P„ ~ Poisson(ia^2") and E denotes expectation. Now note that the expres- 
sion (^/i/2"+^)/sinh(^ft,/2"+^) is the characteristic function corresponding to the 
random variable Logistic(/i/2"+^7r) ^ (ft,/2"+^7r) • Logistic(l). The result now fol- 
lows by interpreting the product of two characteristic functions as the sum of two 
corresponding independent random variables. D 

Hence our proposed simulation method for the Levy area is as follows. Truncate 
the logistic expansion in Theorem 12.11 to include the terms n = 0, 1, . . . , A^. Hence 
consider the approximation AN{h) to the conditioned Levy area A{h) given by 

N 

AnW 



^ n=0 fc=l ^ 



2tt 

where for n = 0, 1, . . . ,iV: the P„ ^ Poisson(-ia^2") are independent Poisson ran- 
dom variables, for k = 1,2,..., F„, X — log(t//(l — U) and Xn,k = log(C/n,fe/(l — 
Un.k) with U,Un,k ^ Unif([0,l]) independent identically distributed uniform ran- 
dom variables, and Z '^ N(0,1) is a standard Normal random variable. We now 
consider how to simulate the tail sum. 



Levy area logistic expansion and simulation 



Theorem 2.2. The mean- square error of the Logistic expansion approximation 
AN{h) is exactly given by 

E\A{h)-A^ihf = j^^('^' 



If we supplement Ai^{h) by including the term (a/v 3 • 2^+^) {h/2) Z to simulate the 
tail sum, where Z ~ A/(0, 1), then the mean-square error in the Logistic expansion 
approximation with tail simulation is 



A{h) - [ AN{h) + , ° J Z 



E 
Proof. Wc directly compute 

/I ^^ -1 ^ Tl 

E\A{h)-A^ihf^E(- y: ^e^^ 



^ 8 1 //i^ ^ 



15 22^+2 V 2 



P„ ^ 2 

27r ^^ 2"^^"'"'' 

^ n=W+l fc=l 

1^) E iEp{^"-^}-nE^-^ 

n=N+l 1=0 ^fe=l 






n=Af+l i=0 



E ^-^{^4 

=Af+l 
oo -. 

2^ 02n 



3V2/ ^^ 22 



3V2/ ^ 22 



n=Af+l 



2on-l 
a 2 



3 \2J 2^+1' 

In order to obtain the error bound for the Normal tail sum approximation, we note 
first that the tail sum has class G distribution. Recall that an infinitely divisible 
random variable has class G distribution, if its characteristic function has the form 
0(^) = exp(-*(^2))^ ^j^gj.g ^(0) = 0, and (-l)"-i*(")(^) > for all n, see Rydcn 
& Wiktorsson (page 163). Using that 

A{h)-A,ih) = ± Y. ^E^-'^^I^Ef^^"'^' 

n=N+l k=l n=0 fc=l 

where h = /i/2^+^, and where Qn has a Poisson distribution with parameter ^a'^2'^ 
and where a^ = a'^2^~^^, it follows that the tail sum A{h) — Ajv(/i) has the charac- 
teristic function (see Levy 1951) 



/ 00 

,^jv(C) =exp('-ia2(i/i^coth(i/i^)- l)) =cxp| -ia^^ 



e 



\ {nh/2TTY + e, 
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Hence A{h) — A]\f{h) has class G distribution — see also Proposition 5 in Rydcn & 
Wiktorsson (2001). We can now proceed as in the proof of Theorem 7 in Ryden 
& Wiktorsson. The tail sum can be represented as a product of a standard Nor- 
mal random variable Z and the square root of an independent positive, infinitely 
divisible variable Yn, i.e. A{h) — A]s[{h) = Z^JY-^. If u\ denotes the variance of 
Aiji) — Ajsjiji)^ then the mean-square error when including the Normal tail approx- 
imation is given by 



E|^(/i)-ylw(M 




Let g denote the Laplace transform of Yat. Then g{z) = (^j^{\f2z\ and the variance 
of Y]\s is given by (logg)"(0). If t[z) :— — (^^/2z coth{^/2z) — l) then we see that 

logg(z) = -a^fz^) and (log5)"(0) = ~a' (^) £"(0). 



Thus we see that 



2^1 ^2^^^ "'"-' 8 1 ^ ^^ 



E\Aih)-A^{h)-a^Z\ ^^^\^) ^"(0)^15 22^+2 ^2 
giving the required result. D 

Note that our expression for the mean-square error of our Logistic approximation 
Aff{h) is exact. Our proposed Logistic approximation to A(h), including simulation 
of the tail sum is thus AN{h) + {a/VS -2^+1) {h/2) Z, where Z - N(0, 1). At each 
order n in the Logistic approximation we must on average sample E{P„} — ^a^2" 
uniform random variables. Hence the complexity of the Logistic approximation 
A]\[{h) is £~^, while with tail simulation it is e^^'^. 

3. Normal approximation 

We have seen that when we measure the complexity of a simulation method in 
terms of the number of uniform random variables required to achieve a given ac- 
curacy, the Logistic, Laplace (Ryden-Wiktorsson) and Levy expansion simulation 
methods scale in the same way. Here we attempt to exploit the form of the Logistic 
expansion simulation method to significantly improve its complexity competitive- 
ness. For each n ^ in the Logistic expansion we must generate and accumulate P„ 
uniform independent identically distributed random variables where P„ is a Pois- 
son random variable with mean E{P„} = ^a^2". Thus on average the number of 
uniform random variables we must generate and accumulate grows exponentially 
with n. One possibility to achieve enormous reductions in computational effort for 
large n, is to proceed as follows. Suppose we generate P„ and realize P„ ^ P for 
a predetermined P ^ 1 to be defined presently. Then in this case we approximate 
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the Logistic random variable accumulant by a suitable Normal random variable. 
More precisely, for each n ^ for which P„ ^ P we approximate 



2_^^n,k ~ V-Pn/37rZ„, 



fe=l 



where the Z„ ^ N(0, 1) arc independent standard Normal random variables. Using 
the Berry-Essecn Theorem we can derive a uniform estimate of the error in the 
distribution function due to this replacement. 

Theorem 3.1. Let Y/v denote the random, variable for some given P G N; 

\ ,1=0 ^fc=l 

and let Fy denote its distribution function. Let Fa denote the distribution function 
of AnIH). Then there exists a constant c such that for all a; G R, 

\Fa{x) ~ Fy{x)\ ^c- 



P 

Proof. The assertion will follow from the Berry-Esseen Theorem. Indeed, let P„ 
denote the distribution function of (1/2") X]fc=i ^n,k ■ 1{p„<p} + \JPnl'^ ■ t^ ■ Zn ■ 
^{Pn^P} ^^^ /" it^ probability density function, and let Gn denote the distribution 
function of (1/2") X]fe=i -^n,k and gn its probability density function. Thus the 
distribution functions Fy and Fa are given as follows (here '•' denotes the usual 
convolution between probability density functions corresponding to the sum of two 
random variables): 

/27Tx/h 
ifx*fi*---*fN){y)dy, 
-OO 

and 

/27Tx/h 
{fx*gi*---*gN){y)dy, 
-OO 

where fx denotes the probability density function of the standard Logistic random 
variable X . The distribution function Fn is given by 



Fnix) 



^ P{P„ = /} • P i ^ X,„k ■ l{i<p} + VU^nZ,, ■ l{,^p} «; x2" I 

1=0 ^k=l ^ 

P-1 f I N OO 

^P{P„ = 1} • P| ^X„,fc < a;2" U^P{P„ = /} • P{ VV3^Z„ ^ a; 2"}, 



and the distribution function Gn by 



OO f I X 

Gn{x) =J2^{Pn = 1} -F \j2^n,k ^ X ■2"\. 
1=0 ^k=l ^ 
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By applying the Central Limit Theorem and the Bcrry-Esseen Theorem to the sum 
of logistic random variables X„.fc on the right, we know that there exists a constant 
c, independent of n, such that for all I and for all x 



\F„{x) - Gnix)\ 

i=p ^ '-fc^i '' 

oo 
l=P 

Since | Fa (x) — Fy (x) | is given by 

2Trx/h 

(/x * (.91 - /i) * 52 * ■ ■ • * ffAf + • • ■ + /x * /i * /2 * ■ • • * (ffjv - In)) (y) dy 
it follows that 

N N oo 



=0 n=Ql=P 

where we have used the following inequality for the convolution product 



\Fa{x) - Fy(x)| < ^ ||G„ - F„||oo ^ C ^ J2^{Pn = 1} / Vl < C ^- 



{f*9){y)dy 



^ \\G\\ 



for any probability density functions / and g, where G{x) ~ J_ g{y) Ay. D 

The Logistic Normal approximation simulation method we propose is the ran- 
dom variable Y^ih) defined in Theorem 13. II We can further augment this method 
by also simulating the tail, i.e. by simulating iV(^) + (a/\/3 • 2^+^) {h/2) Z , where 
Z ^ N(0, 1). Theorem 13.11 demonstrates that the error in the distribution function 
corresponding to A{h) as a result of making the Normal replacement Yis[{h), has an 
upper bound that scales like {N + \)/\/P. In practice, as we see from our simula- 
tions in Figure [21 this upper bound is somewhat pessimistic. On the other hand, if 
we were to be overly optimistic and supposed that this Normal replacement (of the 
large sum of Logistic random variables) did not impact the accuracy of the simula- 
tion, then we would conclude that above a certain computational effort threshold, 
the complexity of the Logistic Normal approximation simulation method would 
be log £ — we only need one uniform random variable for each Normal sample and 
the terms at each order scale like 1/2". In reality the complexity lies somewhere 
between £~^ (or e~^'^ when we include the tail approximation) and loge. 

4. Exponential product approximation 

A Logistic random variable can be decomposed as the difference of the logarithms 
of two independent Exponential random variables. We exploit this result to derive 
a highly effective sampling method for a large sum of logistic random variables. 
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Indeed, the sum of P independent standard Logistic random variables can be rep- 
resented by 

p p p 

^ Logistic,.(l) ^ logYn Exp,(l)) - logffl Exp^.(l)Y 

Hence we are left with the question of whether we can efficiently sample from the 
probability density function associated with the logarithm of the product of P stan- 
dard Exponential random variables (with unit means). In fact the density function 
for the logarithm of a product of P Exponential random variables can be computed 
via the the inverse Mellin transform as follows (see for example Nadarajah 2009; 
p. 654). We briefly recall some basic results on products of random variables and 
the Mellin transform. If two random variables X and Y have densities jx and /y , 
respectively, then their product XY has density given by the product convolution 



f fx(-)fYiy)-dy. 
Jr \yj y 



The Mellin transform of a function has many properties analogous to that for the 
Fourier and Laplace transforms — see Flajolet, Gourdon and Dumas (1995) for more 
details. Importantly the Mellin transform of the product convolution of two func- 
tions is given by the real product of their Mellin transforms. Hence since the density 
of the product of P Exponential random variables is the P-fold product convolution 
of the corresponding densities, its Mellin transform is the P-fold real product of 
the Mellin transform of the standard Exponential random variable density which is 
the Gamma function. The probability density function ip = ip{x) corresponding to 
product of P standard Exponential random variables is thus given by the inverse 
Mellin transform of the Pth power of the Gamma function defined for all a; > by 

-1 />c+ioo 

^(x) = — / x-^{T{z)fdz, 

where c > is any real constant. The probability density function (jjp — (j)p{x) 
corresponding to the logarithm of such a product of Exponential random variables 
is given by (pp{x) — <p(e^) e^. In other words we have 

/•c+ioo 



-1 /"C-^lOO 

P(x) = — / e-(i-^)(r(z))^dz. 



We can explicitly compute the contour integral above using the theory of residues 
from complex analysis to derive a series representation for the density (j)p; see 
[Appendix A[ Unfortunately the coefficients of this series grow so fast with P, that 
its radius of convergence and range of practical evaluation makes it redundant for 
values of P greater than 10. We can also derive the leading order asymptotic form 
for (j)p for large x. However again unfortunately, for large P, the range of large 
values of x for which this asymptotic form is useful, becomes severely restricted. In 
practical terms, the most useful representation is the leading order asymptotic form 
for (f)p for large P, uniform in x, derived using the method of steepest descents. 
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Theorem 4.1 (Logarithm of Product Exponentials density). For large P, uni- 
formly in X, the probability density function (j)p — (f>p{x) for logarithm of the product 
o/ P e N Exponential random variables, has the leading order asymptotic form 



(l)p{x) 



^2t:PiI)'{iI)~^{x/P)) 



where ip = F'/F is the Digamma function. The next order term is asymptotically 
smaller than a factor \/P times a functional form similar to the numerator for (j)p 
above (rather than the factor 1/ yP shown for 4>p)- 

Proof. To determine the asymptotic behaviour for (pp for large P, we apply the 
method of steepest descents. Recall that (t>p = 4>p{x) is given for any real c > by 

0p(a;) = ^ / e-(i-) (F(z))^dz = — / e^((-^/^)^+^°sr(-)) dz. 

27" Jc-ioo 27ri J^_i^ 

Let /i(z) := {—x/P)z + logF(z) denote the exponent function in the integrand. 
Note that we retain the ratio {—x/P) as we wish to include the possibility that 
(—x/P) might be order one or even asymptotically large. Since P and x are real, 
we note that h — h{z) has a unique saddle point given by z = z* := ip~^{x/P) > 
where ip = F'/F is the Digamma function. Since the argument in the above contour 
integral is analytic in the right-half complex plane, we can deform the contour 
of integration to one of constant phase that passes through z = z* on the real 
axis. Then we apply the usual Laplace method arguments. We shrink the range 
of the contour of integration to a small interval around z = z^, replace h{z) by 
its quadratic Taylor expansion, then extend the range of integration out again and 
utilize standard area estimates for Gaussian functions. This standard procedure 
reveals that 



6' 



9_ /„P/i(z.) 



^^ ' 27ri y Ph"{z^,) V P 

where h"{z^,) denotes the second derivative of h along the constant phase contour 
(orthogonal to the real axis) at z = z*. Substituting the form for h — h{z) and that 
z* = tlj~^{x/P) gives the large P result of the theorem. D 

As with the Logistic Normal approximation, the idea here is to more efRciently 
sample the large sums of Logistic random variables, on average E{P„} = ^a^2", 
at each order; especially at higher orders. Hence when the Poisson sample P„ ^ P 
at order n, for some predetermined P ^ 1, we will instead draw samples from the 
distribution function ^p^ for the logarithm of P„ standard Exponential random 
variables — corresponding to the density function cjjp^ . The only practical form for 
(j)p^ for sampling when P„ > 10 is the asymptotic form given in Theorem 14. II How- 
ever, unfortunately an analytic form for distribution function $p^ is not available. 
To overcome this we proceed as follows. We tabulate the corresponding inverse 
distribution function ^^^i for £ — 2,3,.... This needs to be done only once and 
as accurately as we possibly can. Then given any P„ ^ Poisson (ia^2") such that 
P„ ^ P, with say P = 10^, we can decompose P„ = ai -I- 02-10^ -I- as •10'^ -I- • • • . Here 
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ai,a2, . . . are simply the decimal digits of P„. In the Exponential product approx- 
imation method we propose then, for each n, when P„ ^ P we sample ai Logistic 
random variables Xi , . . . , Xa-^ directly. Then using interpolation with the table 
for $~p2 we sample 02 pairs of random variables C/2,1, . . . , C/2,a2 a-nd V2.1, ■ ■ ■ ,V2,a2 
from $102. We further sample 03 pairs of random variables C/3,1, • • • , C/3.a2 and 
V3.1, . . . , V3_a2 using interpolation from our table for ^^qs, and so forth. Hence for 
Pn ^ P, we replace the Logistic random variable accumulant as follows 

Y^ XnM ~ X! ^J + X! X! ^''^0 ~ "^Lj- 
k=l j = l £^2 j = l 

In practice we restricted ^ ^ 5. The only error in this replacement is in the asymp- 
totic approximation of the cumulative distribution function, the accuracy of the 
tables and the interpolation. Above the threshold P, the number of uniform ran- 
dom we require at each order is two times ai -|- a2 -I- 03 -|- . . .. Hence, since at each 
order the Logistic accumulant is scaled by the factor 2", the complexity of our 
Exponential product approximation method is in principle loge. We can of course 
augment this method by simulating the tail sum as before. 

5. Simulation details 

We provide here details of our implementation underlying Figure [5] and at the same 
time explain more carefully what we observe therein. Our Levy and basic Logis- 
tic simulation are straightforward truncations of the corresponding representations 
given in the introduction and Theorem 12.11 Note that in all cases we performed 
L — hx 10® simulations — as large as was practically possible. Monte-Carlo simula- 
tion variance is of order 1/L. From our experience and as can be seen in Figure [5J 
there is significant Monte-Carlo simulation noise roughly below 10~^, which we 
should take into account when intepreting the mean-square errors therein. Note 
that we measure the error as the absolute value of the sample variance minus the 
true variance (1 + a^) /i^/12, where with random Wiener increments /SW^{h) and 
AW^{h) for each sample, we set a^ = ¥.{{l^W^{h)f + {l^W'^{h)f] = 2h. We see 
in Figure [21 that the Levy method provides Levy area samples for low accuracies 
and low computational effort. The slope of the error versus the actual computa- 
tional effort for the Levy method is —1 as it should be (the variance of the tail 
asymptotically decays like 1/n). If we also simulate the tail for the Levy method, 
then the mean-square error should decay like 1/n^; see Kloeden, Platen & Wright 
(1992) or Wiktorsson (2001). Note that we do not directly observe this improved 
accuracy/convergence in Figure [H as we compute the sample variance for the Levy 
method with independent tail approximations. Thus we expect to see the error of 
the Levy method with tail approximation to be similar to the error of the Levy 
method itself, shifted by the variance of the tail approximation — an independent 
computation confirms this. Indeed such a shift is observed for all the methods we 
implemented when we also simulate the tail, except that the shift is much more 
dramatic for all the Logistic expansion methods. In the case of the basic Logistic 
method, we can only generate Levy area samples if we're willing to invest compu- 
tational effort above a (in principle quantifiable) threshold. In the basic Logistic 
approximation, the order scales exponentially — the terms at each order scale like 



14 Simon J. A. Malham and Anke Wiese 



2^" but we must invest on average a factor E{P„} = -^0^2" amount of effort 
at each order as well. In Figure [2l low order basic Logistic approximations scale 
favourably as it's relatively cheap to sum a medium-sized number ia^2" of Logistic 
random variables (whilst the terms scale like 2~" at each order). However eventu- 
ally, as we can see in Figure [51 we reach the asymptotic behaviour we expect with 
a slope of —1. When we also simulate the tail, we get a dramatic improvement 
in accuracy, so much so that the results lie in the regime where we know there is 
significant Monte-Carlo simulation noise. Note that we did not include in Figure [5] 
simulations of the Ryden-Wiktorsson method. The computational effort required 
to achieve the corresponding accuracies we have indicated in Figure [5] meant that 
the Ryden-Wiktorsson method, even including simulating the tail, was not com- 
petitive. Note also that for all our Logistic expansion based methods, we used the 
standard algorithm given in Knuth (1998) to simulate the Poisson random variable 
Pn when the mean ia^2" ^ 10^ and the Normal approximation otherwise. 

In Figure m we also indicate the performance of our Normal approximation and 
Exponential product approximation methods. In the case of the Normal approxi- 
mation, at each order, when the Poisson sample Pn > 10^ we invoked the Central 
Limit Theorem and replaced the Logistic sum by the suitably matched Normal 
random variable as indicated in fJJ] This effect is strikingly indicated by the sharp 
downturn in the Normal approximation plot in Figure [51 We also observe that the 
accuracy of the method is robustly uncompromised by this Normal replacement. 
Note that including a tail approximation generates a much more accurate sampling 
method. As we indicated above, an independent computation demonstrates that 
the error shown in Figure [51 for the Normal approximation with tail is the same 
as the Normal approximation shifted by the variance of the tail approximation. 
Note that the Normal approximation with tail carries us into the regime where the 
Monte-Carlo simulation noise is significant. 

To implement the Exponential product approximation some preparation is re- 
quired. We constructed tables for the inverse distribution functions $p for the 
logarithm of the product of P Exponential random variables for P = 10^, lO'^, 10^ 
and 10^ as follows. Using the large P asymptotic approximation for ipp we derived 
in Theorem 14.11 we constructed ^p by integrating from the left and right ends, 
respectively, as follows (we justify this presently): 

/X I'OO 

(l)p{y)dy and ^p{x) = 1~ <j)p{y)Ay. 

-oo J X 

We used the trapezium rule with very small stepsizes to approximate these two 
expressions and the results are given in Figure [31 In the left set of panels of the 
Figure we plot the distribution functions <i?p constructed in this way. The blue 
curve corresponds to integrating from the left and the green to integrating from the 
right. An empirically computed cumulative distribution function is shown in red. 
The accuracy of the integrated asymptotic form is immediately apparent; as is the 
reason for independently integrating from either end to minimize the accumulated 
quadrature error. For the case P — 10^, the magenta curve shown corresponds to 
constructing the distribution function by a quadrature approximation of the contour 
integral form for <^p — $p(a;) for each x. For small x this approximation breaks 
down and for larger values of P it becomes impractical. In the right set of panels 
in Figure [31 we show the difference of the left and right integral approximations 
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-5.9 -5.85 -5.8 -5.75 -5.7 -5.65 

" xlO* 



Figure 3. The panels on the left show the asymptotic approximations (acdf: from the left, 
blue, from the right, green) and the empirical cumulative distribution function (ecdf: red; 
for which we used a 'stairs' plot) for the logarithm of the product of P standard Exponen- 
tial random variables. Note we have used a log-of-minus-log scale. The panels on the right 
give an indication of the error of the the asymptotic approximations compared to the em- 
pirical cumulative distribution function (with the same colour attribution). The magenta 
curve in the top two panels represents an approximation to the cumulative distribution 
function constructed by directly numerically computing the contour integral for (/)p. 
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to the empirically computed cum.ulative distribution function which indicate errors 
less than 10^'^. We spliced together the left and right approximations, numerically 
inverted the result, and then used linear interpolation to construct the inverse 
^p^ ~ ^p^{x) cumulative distribution functions using 10^ + 1 points uniformly 
spaced on [0, 1]. In other words, the tables each contain 10^ + 1 values for ^p^ at 
Xm = rn/10^ ioi m = 0,1,2, ... , 10^. In the tables corresponding to P = 10'^ and 
above, we crudely set $p^(l) — 1 whilst we set $"^2(1) = 10. Given the range 
of values of $p^ as x — >■ 1 we can infer from the left set of panels in Figure [3] 
this replacement of +00 seems reasonable and in practice from what we observe 
in Figure [5] this did not affect accuracy. The tables thus constructed in this way, 
allow us to rapidly sample from <I>p for any of the values of P we've indicated, as 
follows. We generate a uniform random variable U. The integer part oi U x 10^ 
immediately indicates the position in the table of interest to us. We use linear 
interpolation between that position and the next to generate a sample from $p. 
Thus, with these tables and linear interpolation method for sampling from them, 
we implemented the Exponential product approximation as indicated at the end 
of 21 We can of course also include a tail approximation. The results shown in 
Figure [5] are just as dramatic as for the Normal approximation, making both these 
methods with tail simulation the ones of choice for high accuracy simulations. 

To conclude we remark that the improved sampling methods for the Levy area 
we have introduced have important implications for our theoretical work on the 
algebraic structure of stochastic differential equations and concomitant universally 
accurate integrators (see Malham & Wiese 2009 and Ebrahimi-Fard et. al. 2011) as 
well as for more practical work on the pricing of financial derivatives (see Malham 
& Wiese 2011). 

Appendix A. Series representation 

Applying the theory residues to the contour integral form for the probability density 
function (jip, corresponding to the logarithm of the product of P standard Expo- 
nential random variables, we can explicitly derive a series representation for it. 
For convenience, we introduce the discrete convolution product. For two sequences 
a — {ao, oi, 02, . . .} and b = {b^, 61, &2, . • ■} this is the sequence whose nth term is 
given by (a * 6)„ := a„&o + a?i-i^i + • • • + oofon, for n = 0, 1, 2, . . .. If a and b are the 
coefficients of the corresponding powers in a power series, say oq + aix + a2X^ + • • • 
and &o + bix + b2X^ + • • • , then the coefficients of the product of these two power 
series are a * 6. 

Theorem 5.1 (Series representation). For the probability density function <j)p — 
(j}p {x) for the logarithm of the product of P Cz N Exponential random variables we 
have the following explicit series representation 

p-i 

0p(a;) = ^a;''^^afc(n)e("+l)^ 

fe=0 n>Q 



where 






'^fc^") - (^UPU ■ (9 * ^W)p-l-fc- 
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Re(z) 



Figure 4. The contour Cb, with sections Cji, C\ and Cj{ used in the proof of Theorem 14.1 



Here ifT^^'{l) is the kth derivative of the Gamma function evaluated at one, and 
if C^ denotes the usual Binomial coefficient, P choose k, then 



,:=({r('=)(l)/fc!},J 



*p 



and 



K-)-{CfUo*{^f/2'} 



k>0 



{C^/n'] 



k>0' 



Proof To derive the power series representation for the probabihty density function 
(j>p, we compute the inverse Mellin transform ip of F^. Here ip — (p{x) corresponds to 
the probabihty density function of the product of P standard Exponential random 
variables. At the very end we use that (i>p{x) = ^p{e^) e^. 

To analytically compute the inverse Mellin transform of F^ we use the theorem 
of residues for the contour Coo = limfl^cxjC/? shown in Figure S] we found Gajjar 
(2010, p. 111-135) very useful here. Note that we have 

-1 /•c+ioo -I /• 

= ^res(x-"(F(z))^: z = - 



Here we have used that the contributions to the contour integral from the com- 
ponents of the contour C]^, C^ and C^ are zero. Our goal now is to determine the 

residue of a;~^(F(z)) at each of its poles z — —n for n = 0, 1, 2, Note these 

are poles of order P due to r{z) have simple poles at these points. Our goal now 
is thus to derive the Laurent series in powers of {z + n) for x~^ (^(z)) centred at 
each of these poles, and to determine the coefficients of (z + n)^^ corresponding to 
the residues. For convenience we set C, := z + n. The Taylor series for x~^ centred 
at z = — n is given by 



• exp(— Clogx) — x' 



E 



(-logx)*^ .fc 



A:! 



C" 



E 

fcSO 



4■(x)c^ 
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where we set £k{x) ~ (— logx)'^/fc! for all k ^ 0. Then using the Taylor series for 
r(l + z) centred at z = 0, we have 



r(z) = z-\z + 1)-^ • • • (z + n)-^T{l + n + z) 

■} 

~k\ 



(c-nr(c-n + i)-^---r^-E^c^ 



fc>0 



Taking the P power of this expression, and using the definition for the coefficients 
{qk\ and {rk{n)} given in the theorem, we get 



(r(.))" = (c - n)-'{c - (n - 1))-- ... (c - i)-r- • (e ^C^ 



p 



(i-c/n)""(i-c/(n-i))""---(i-c)"V^-E*c'= 



(-1) 



k>0 



Pn 



C^-J2{r{n)*qU''- 



{nl)P 

Thus we are now in a position to compute the Laurent series for x^'- (r(z)) centred 
at C = corresponding to z = —n. Directly computing, we have 

l_\\Pn , , 

= 7^--"-E(^(-)*(-H*'^)),-C'^"'^' 

where f(x) = {i!fc(a;)}fc^o- The residue of x^^(r(z)) at C = is generated by the 
term fc = P — 1, i.e. we have 

res(x-^ (rlz))"^: z - -nj - ~uh^ ' ^" ' (^^^^ * ^("^ * '^)p-r 

Using the explicit sum for the first discrete convolution on the right-hand side, 
summing over all the residues z = 0,— 1,— 2,..., and using that 4>p{x) — (^(e^) e^, 
generates the form for the density function i/ip given in the theorem. Note that 
r(0) = {1,0,0,...} and so {£{x) * r{0) * q) p_^ = £o{x)qp-i + h{x)qp-2 + ■ ■ ■ + 
£p-i{x)qQ. D 

This explicit series representation is useful and accurate for relatively small 
values of P. In particular (f>2 matches the known modified Bessel function form 
2Ko(2e^^^)e^ for the probability density function when P = 2. Further, integration 
by parts soon generates an analytic series expansion for the corresponding distri- 
bution function. Unfortunately however, the coefficients q grow factorially with P 
making the series impractical for P greater than 10. 
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As mentioned in the main text, we can also derive the leading order asymptotic 
form for (^p for large x. Consider the product convolution form for density of the 
product of P standard Exponential random variables given by 

, , roo roo . X ^ exp(-(y2 + --- + ;/p)) ^ , 

tp{x) ^ •••/ exp dy2...dyp. 

Jo Jo \ y2---ypj V2---yp 

After using the substitution yi — UiX^'^, we can then apply Laplace's method. 
Noting that the exponent function has a global maximum at 2/2 = • • • = yp = 1 in 
the range of integration, then up to exponentially small errors, we have that 

^ ^ exp(^Px^/P) r°<' r°<' , . . , \- 1 \ 

'^E(^) ^(p_l)/2P • J^ ■ J^ exp( - > ^ T, - > ^ ^nTj ) dn . . . dTp_l. 

The matrix generating the quadratic form in the exponent in the integrand has one 
eigenvalue equal to P — 1 with the remaining eigenvalues all equal to one. A linear 
change of coordinates corresponding to the orthonormal similarity transformation 
of the matrix generating the quadratic form, and using that (j)p{x) = (p{e^)e^ 
generates the asymptotic form 

, , p-i 

(l>p{x) ~ ^ ^j_' • exp(-Pexp(x/P) + {P + l)x/2P). 
vP 



p-l 


p-i 
- > ] w, 



For large P though, the range of large values of x for which this approximation is 
useful, is severely restricted. 

Lastly, we remark that an analytical form for the probability density function 
corresponding to a sum of Logistic random variables is given in George and Mud- 
holkar (1983). We have not pursued this approach, though one could in principle 
investigate whether a computationally feasible method to sample a sum of Logistic 
random variables could be derived from this form. 
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